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Abstract 

New resonance steps are found in the experimental current-voltage charac- 
teristics of long, discrete, one-dimensional Josephson junction arrays with 
open boundaries and in an external magnetic field. The junctions are under- 
damped, connected in parallel, and DC biased. Numerical simulations based 
on the discrete sine-Gordon model are carried out, and show that the solutions 
on the steps are periodic trains of fluxons, phase-locked by a finite amplitude 
radiation. Power spectra of the voltages consist of a small number of harmonic 
peaks, which may be exploited for possible oscillator applications. The steps 
form a family that can be numbered by the harmonic content of the radiation, 
the first member corresponding to the Eck step. Discreteness of the arrays 
is shown to be essential for appearance of the higher order steps. We use 
a multi-mode extension of the harmonic balance analysis, and estimate the 
resonance frequencies, the AC voltage amplitudes, and the theoretical limit 
on the output power on the first two steps. 
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I. INTRODUCTION 



Josephson junction systems have a natural application as millimeter- 
and submillimeter-wave oscillators. To facilitate their use, however, ongo- 
ing research must address several issues. To produce oscillators with narrow 
linewidths and high power, we must study the harmonic content and oscil- 
lation amplitudes of Josephson sources. Furthermore, the conditions under 
which many Josephson oscillators can be phase-locked to produce higher out- 
put power continue to challenge researchers. 

AC power has been measured from both continuous long Josephson junc- 
tions and discrete arrays of short junctions. In underdamped Josephson sys- 
tems, the AC oscillation amplitudes are expected to be largest at certain 
resonant frequencies. Much analytical work has been done to predict the fre- 
quency of these resonances. The success of these analyses can be assessed by 
studying the DC current-voltage (I—V) characteristics. Steps appear in the 
I—V when increases of the current bias over a certain range do not produce 
increases in the DC voltage. Instead, the input power drives large-amplitude 
AC oscillations. The frequency of these oscillations is related to the DC volt- 
age by the Josephson relation, u> = 2irV/& , where <E> is the flux quantum, 
while the spatial wavenumber is determined by the applied magnetic field, 
k = 27r/, where / is the applied flux per unit cell normalized to The os- 
cillations can be approximated as the normal modes to the linearized system. 
This approach has been used successfully to predict resonance frequencies, 
or steps in the I—V characteristics of many different geometries of Josephson 
systems. The dependence of the step voltage on applied magnetic field can be 
included, thus giving a dispersion relation for the associated oscillations. In 
the discrete case, the dispersion relation is nonlinear, and higher modes waves 
are expected to produce distinct steps in the I-V [Q^). One of the limita- 
tions of these linear analyses is that the predicted oscillations are necessarily 
single-harmonic, since the normal modes of the linearized Josephson system 
are simply Fourier modes. In addition, the absence of a driving force in the 
linearized system precludes any calculation of the oscillation amplitudes. 

In order to predict oscillation amplitudes and the distribution of power 
among excited modes, the nonlinearity in the system must be more carefully 
treated. Perturbation techniques (^,^| clarify the role of the nonlinearity in 
driving the resonance and have been used to predict step voltages in Josephson 
systems. Combined with an appropriate ansatz, they also provide expressions 
for the oscillation amplitude [f^-Q]. However, analytic expressions for the AC 
amplitudes so far include only the first harmonic and are implicit ||. 

In this paper we use experiments, simulations, and non-linear analysis 
to investigate the properties of discrete, planar arrays of Josephson junctions 
connected in parallel. Both experiments and simulations yield several resonant 
steps in the I—V. Simulations indicate that on these steps, a traveling wave 
pattern dominates the oscillations of junctions in the array. In contrast to the 
linear picture, these resonances correspond to excitations of more than one 
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FIG. 1. Schematic of an open-ended parallel array. A uniform current lb is applied at each of 
the upper nodes and extracted at the lower nodes. 



mode, though the harmonic content is still limited. Similar states have already 
been reported in discrete rings ||, indicating that boundary conditions play a 
minor role. In inductively coupled arrays with open ends, the second harmonic 
resonance was also observed experimentally Q. We use a two-mode extension 
of the harmonic balance method to predict the resonant frequencies as well as 
the mode amplitudes at resonance. Our formulas are analytic and include the 
dependence on magnetic field. With a matched load condition, a theoretical 
upper limit on the available output power of the underdamped Josephson 
oscillator is calculated. 



II. MEASUREMENTS OF ARRAYS 

We have measured single-row arrays of N = 54 junctions connected in 
parallel. Figure [l] shows a schematic of our device. A bias current lb is 
applied at each upper node of the array and extracted from each bottom 
node as shown. We use resistors to distribute the bias current as evenly as 
possible. The voltage across the array is measured at an edge. The array 
is placed above a superconducting ground plane, and a separate control line 
(not shown) is used to apply a magnetic field. We will discuss the applied field 
in terms of the frustration, f (the flux applied to a single loop of the array, 
normalized to the flux quantum). Because the system is discrete, we expect 



its properties to be periodic in / with period / = 1 [10|. In this experiment, 
the applied flux is proportional to the control current. 

Samples were fabricated using a Nb trilayer process filf . The junctions 
are 3 x 3 fim 2 with a critical current density of j c (T = 0) = 1270 A/cm 2 . 
Device parameters have been determined using the diagnostic procedures de- 



scribed by van der Zant et al. |12j. For our samples at 7.2 K, the normal- 
state resisitance R n = 16.6 0, the self-inductance of a loop L s = 6.4 pH, the 
nearest-neighbor inductive coupling = 0.11L S , the junction capacitance 
C = 340 fF, and the Josephson inductance Lj = $ /(27r/ c ) = 5.9 pH. We use 
the normal-state resistance to calculate the Stewart-McCumber parameter, 
Pc = 16. The sub-gap resistance is not as well defined at these high tem- 
peratures but approaches the value of R n . For the discretness parameter, we 
calculate A 2 = Lj/L s = 0.92. 
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FIG. 2. Current vs. voltage of a 54-junction array on a ground plane at three values of 
/ = 0.2,0.3,0.5. The temperature is 7.2 K so that I c R n = 0.93 mV. The three steps indicated 
by the arrows are labeled by m values corresponding the number of dominant harmonics in the 
mode. 



0.5 




0.2 0.4 0.6 0.8 1 



f 

FIG. 3. Step voltages of the 54-junction array vs. the frustration, /, for the modes m = 1 and 
m = 2. The dashed curves are plots of Equation (13). 

Figure || shows the voltage across an array when the current is uniformly 
injected. Three I—V curves are presented, for / = 0.2, 0.3, and 0.5. The 
most prominent feature is the sharp Eck step, which is the steepest part of 
the I—V, just before the switch to the gap voltage (not shown). We label this 
step "m = 1" for the reason given in the next section. In the regime where 
Pc is small (overdamped) or Aj is large (less discrete), this is the only step 
observed. For our underdamped and discrete samples, however, additional 
steps appear below the Eck voltage. We index these steps with m = 2,3, etc. 
In this paper we are mostly concerned with the Eck step m = 1 and, among 
the new steps, the clearest one with m = 2. The study of the second step 
sheds light on the finer steps with m > 2. 

The voltage locations of the peaks experimentally vary with /, as seen in 
Fig. |2[ More systematically, we show the dependence of the first two steps in 
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Fig. ||. The Eck peak voltage is found to be periodic in / with period / = 1 
and to be approximately symmetric with respect to / = 0.5. This is consistent 
with previous observations [Hj. At / = 0.5, the Eck step reaches its highest 
voltage value. For a smaller /, there is a threshold frustration, below which 
the Eck step does not appear. This cut-off f c \, known as the lower critical field 
or frustration, is the minimum applied flux density for vortices to enter an 
array [10]. The value is quite large for our system, f c \ m 2/(7r 2 Aj) = 0.2. The 
voltage location of the second step shows roughly the same /-periodicity and 
symmetry as the Eck step. This second step, however, achieves the maximum 
voltage near / = 0.25, and it disappears near / = 0.5 and for approximately 

/ < fcl- 

The Eck (m = 1) steps are ubiquitous in one-dimensional parallel arrays 
as well as in continuous long junctions. In contrast, the other steps (m > 1) 
do not appear in long continuous junctions, but do appear in discrete arrays 
when Aj is small. In our arrays, we find Aj must be less than unity for 
the m = 2 step to appear. Similar steps have been observed also in highly 
discrete circular arrays ^ and open-ended arrays consisting of two rows that 
are inductively coupled ||. 

In open-ended arrays with a smaller N, Fiske steps ]l3| may be observed in 
a similar part of the I— V below the Eck voltage. We emphasize, however, that 
they are qualitatively different. Fiske resonances can be described as standing 
waves (cavity modes) resulting from boundary reflections. The wavelength of 
the modes are restricted by the boundary geometry, and consequently, the 
resonance voltage locations do not depend strongly on /. At a certain value 
of /, only even or only odd modes are excited. As N becomes large, for a 
given value of damping, these Fiske resonances disappear due to damping of 
the edge reflections. None of these features apply for the Eck step as well as 
the m > 1 steps, which are tunable in /. Thus, the new steps are expected 
to belong to the same family as the Eck step. 



III. SIMULATIONS 

The governing equations which model our arrays are derived by applying 
Kirchhoff's current laws and using the RSJ model for the current through 
a single junction. We normalize the current to I c , the voltage to I c Rn, an d 
time to \/LjC (inverse plasma frequency). The equations are given in terms 
of the gauge-invariant phase differences across the junctions, <f>j, where j = 
1, . . . , N indexes the junction's position. For simplicity we neglect all the cell 
inductances except the self-inductance, which results in the damped, driven, 
discrete sine-Gordon model: 

4>j + T<pj + sin 4>j = I b /I c + A 2 j(cp j+1 - 2(f) j + (1) 

— 1/2 

for j = 1, . . . , N and F = f3 c ■ Our numerical code can include longer range 
inductances as necessary, up to the full inductance matrix p3. However, our 
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FIG. 4. (a) Simulated current-voltage characteristic for a 54-junction array at / = 0.3. The 
three most prominent steps are marked, (b) Voltages corresponding to the top of the steps show 
the tunability of the first two steps with magnetic field. The parameters are the experimental 
values: (3 C = 16 (T = 0.25) and A 2 j = 0.92. 



analysis in the next section uses mainly ([[]). Simulating the same equations 
allows us to make a direct comparison with our analysis. The simpler system 
not only illuminates the essential mechanism responsible for the observed 
steps but also reproduces the measurements reasonably well, as we will show 
in this section. 

With only self-inductance, the boundary conditions are simply 

fo(t) = Mt) - 2?r/ and <p N +i(t) = <j) N {t) + 2vr/ (2) 

for all t, where artificial junctions <f>o and 4>n+i are introduced at the endpoints 
so that (|l|) is valid at j = 1 and N as well [15]. The fourth-order Runge- 



Kutta scheme with a time-step At = 1 is used for integrating the system. 
The instantaneous voltage at junction j is simply proportional to the rate of 
the change of (fij , and is given in our normalization by 

V/IcRn = Tdfa/dt. (3) 

With (0-|3|), current- voltage characteristics are numerically obtained at 
different / values using the parameters Aj and T from our experiments. Fig- 
ure H](a) shows the results for / = 0.3. The curve reproduces the measured 
one in Fig. Q well, and at least three steps (indexed by m) are clear. (There 
is also a small step just below the m = 1 step which we think arises from a 
different mechanism.) In a manner similar to Fig. ^, the /-dependence of the 
step voltage locations is shown in Fig. ||(b). There are slight differences, but 
the main features and voltage values are represented well. 

Simulations allow us to study the solutions in detail. We are especially 
interested in the system dynamics when biased on top of a step. In Fig. |5|(a-c) 
the phase of the junction j = 27 (located in the middle of the array) on the 
m = 1,2,3 step, respectively, is shown as a function of time. The junction 
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FIG. 5. Time evolution of the phase of a junction j = 27 in the center of an array, biased on 
the (a) m = 1, (b) m = 2, and (3) m = 3 steps, respectively. Note the increased frequency content 
as the mode number increases. 



appears to be in a periodic motion. Its phase increases rapidly when a vortex 
(kink) passes by, and it oscillates for the period between passing vortices. In a 
mechanical analog of a junction as an underdamped pendulum a sudden 
overturn of the pendulum is followed by an overshoot, and the pendulum 
"rings" several times until the next vortex passes by. 

It appears in Fig. |5| that the solution on the step m corresponds to m 
such ringings. This can be quantified by studying the harmonic content of 
the voltages @. Fourier spectra of the voltages are shown in Fig. |||(a-c), 
respectively. On the step m, the first m harmonics are dominant, and the 
higher harmonics have rapidly decaying magnitudes. Despite the presence of 
other harmonics, the steps are indexed by the ringing frequency, m. 

Similar plots for the other junctions inside the array appear identical to 
c/>27, except for a certain shift in the time axis. This suggests that the solu- 
tions are well approximated by traveling waves. Near both ends of the array, 
reflections from the ends change this picture. The boundary effects, however, 
decay within 4-5 junctions from the end, and appear to play only a minor 
role in our long arrays. 

Such traveling wave solutions, consisting of a vortex and ringing, have 
been found in circular arrays with the periodic boundary conditions [|,15,17j. 
In these systems, a vortex is trapped in a ring, and as it circulates, it creates an 
oscillatory wake, which phase-locks back to the vortex. It is not surprising that 
a nearly identical situation may arise in a linear array with open boundaries, 
if the array is sufficiently long. Instead of a single circulating vortex, a vortex 
lattice propagates through the array. A junction is swept periodically by a 
vortex in both cases. There is, however, one major difference between the two 
geometries. The magnetic field in the arrayin the array (which controls the 
vortex spacing) can be continuously tuned in the linear geometry while it is 
restricted to multiples of <I>o i n a circular array due to the flux quantization. 
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FIG. 6. Fourier spectra of the voltage corresponding to Fig. |j|(a-c), with (a) step m=l, (b) 
step m=2, and (c) step m=3. Voltages are proportional to the time derivative of the phase in 
Fig. H The average slope of the phase plots corresponds to the DC voltage and is determined by 
the fundamental harmonic. The DC voltage was subtracted before computing the Fourier spectra. 
Note that only m harmonics are dominant for the mth step. 

IV. ANALYSIS 

In this section we seek a more analytical description of the system, and 
estimate DC voltages at each step and the amplitudes of the AC voltage com- 
ponents. The simulations in the previous section suggest that finding traveling 
wave solutions to the governing equation (|l]) is the key to the estimates, but 
this task is not simple in practice. 

The equation (|l|) may be viewed as a variation, in two respects, from the 
integrable sine-Gordon partial differential equation, 



<ptt + sin <p = (p xx 

which possesses traveling kink solutions as exact solutions. Firstly, (|l|) is no 
longer a conservative equation, as it includes the external drive (bias current 
lb) and the loss (T(f)j) terms. The added terms also break the integrability 
and exact solutions are no longer known. A perturbation approach (l^] in 
the conservative limit (lb, T — > 0) has been developed to approximate a kink 
solution. 

Secondly, the spatial derivative in (|j) is discretized in (|l|). In general, non- 
linear wave equations discretized on lattices may exhibit qualitatively different 
solutions from their continuous counterparts, and study of such discreteness- 



(4) 



induced effects is an active current topic in its own right [19|. For sine-Gordon 
systems, in particular, (|l]) was studied by [20,21,1^] without a loss term or a 
drive, and it was found that propagation of a kink introduces a background 
radiation which greatly influences the speed of the kink. As far as we are 
aware, however, there has been no attempt to estimate the amplitudes of the 
induced radiation, especially in the driven case. 
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Strunz and Elmer |l7j recently transformed (|lj) into a system of coupled 
modal equations and pointed out that the superharmonic resonance leads to 
creation of the radiated waves. They expanded a traveling waveform into 
Fourier modes, which may not be the best expansion basis but is a convenient 
one. Since our Fourier spectrum in Fig. ^ shows only a small number of peaks, 
we may truncate most modes and still obtain reasonable predictions. With 
this goal in mind, we review the analysis of [jlT]] in Sec. 4A. The coupling 
terms among modes are truncated and analyzed in Sec. 4B, and we estimate 
the AC voltage amplitudes on the steps. The available power from the array 
is then evaluated in Sec. 4C. 



A. Resonance mechanism and voltages 

We look for traveling wave solutions of the discrete sine-Gordon equation 
([]]) of the form 

<pj(t) = 4>(x) = x + ip(x) 

where 

x = Ut + 2-Kfj 

is the moving coordinate with the wave, and tjj(x + 2ir) = ip{x). The funda- 
mental frequency u is proportional to the DC voltage through the Josephson 
voltage-phase relation (||) while the spatial wavenumber is imposed by the 
external field. If the modulation were absent (ip = 0), then the boundary 
conditions (^) would be satisfied exactly. Since if) is not vanishing, they are 
satisfied only on average, and there should be a correction to (||). We neglect 
this boundary effect and will show that the simplification still leads to good 
estimates of the measurements and simulations. (In a circular array (|2|) is 
replaced by the periodic boundary conditions: (ftj+N(t) = (f>j(t) + 2ttM where 
M is an integer |l5|| . There can be exact solutions of the form (||||) with 
/ = M/N.) 

The periodic function ip can be expanded into Fourier modes 

oo 

^(x)= A me imx 

m=—oo 

with A*_ m = A m . We set Aq = 0, without loss of generality, by shifting the 
origin of time. The phase (j) is an increasing function of x, but the nonlinear 
term sin (p is 27r-periodic and can be expanded as 

oo 

sin </>(*)= 

m=— oo 

with F* m = F m . The coefficients F m can be computed in terms of 
A±i, A±2, . . . from the usual Fourier-Bessel expansions, and therefore pro- 
vide coupling among the modes. By substituting (||]7|,||) into (|l|), we obtain 
a coupled system of modal equations. 
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(S m + imTuj)A m + F m = 



where 



and 



The index m = 1, 2, 
in d) results in 



<5 m = u) 2 m - (muj) 2 , 

0J m = 2Ajsin(m7r/). 
for (^-11). In addition, the balance of the DC terms 

I b /I e = Tu + F . 



The first term on the right hand side is proportional to u, and hence to the 
DC voltage. This term describes the ohmic line in the I—V plane. The second 
term is also a function of u> through the A m 's, and describes the deviation of 
the I—V curve from the ohmic the line. 



The superharmonic resonance [22| may occur in the algebraic system 
(|9,12) and lead to the creation of resonanct steps in the I—V curve. For 
small T, the magnitude of A m would generically become large near 5 m = 0, 
that is, near 



to 



uj m /m, 



m 



1,2,.... 



Then Fq also becomes large, and the I—V deviates substantially from the 
ohmic line near these frequencies (or corresponding voltages). This argument 
neglects the dependence of the coupling terms F m on the complex amplitudes 
A m , but explains our measurements quite well. We plot (13) for m = 1, 2 as 
dashed curves in Figs. ||| and [|(b). The agreement is good without a fitting 
parameter. 

The right hand side of ( |l3| ) is the phase velocity of the rath mode on the 
discrete lattice. The left hand side can be viewed as the vortex velocity. Then, 
the relation can be viewed as the phase-locking condition of the vortex velocity 
and the phase velocity of the mth mode Jl|,|l5| . In the absence of the coupling 
terms F m , a resonance occurs when the vortex lattice moves at the same 
velocity as one of the modes. We note that these resonance frequencies to are 
distinct for different m because of the sinusoidal dependence of the dispersion 
relation ([ll]) on m. If the second order difference in ([!]) were replaced by 
the second spatial derivative as in (|j), the dispersion relation would depend 
linearly on m, and the phase velocities of all the modes would be identical. 
Consequently, there would be only one vortex velocity that can excite modes, 
and we would observe only one resonance voltage. This explains why only the 
Eck step (among this type of step) is observed in a long continuous junction. 
In this sense, our observation of the m > 1 steps is due to the discreteness of 
the system. 
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We also note that, in some parameter regime, the resonating modes may 
be sub-harmonics of the fundamental frequency. Sub-harmonic resonances of 
order n can be sought by generalizing ([<]) to 

oo 
m=—oo 

This is still a traveling wave but has an n times longer period. In a similar 
manner to the above discussion, one may expect resonances when 

= w m,n - (moj/n) 2 i.e. u = nuj m , n /m 

where 

u) m , n = 2Aj sin(2m7r//n). 

The resonances with m = 1 and n = 1,2,3 have been observed by Caputo et 
al. @. 

B. Amplitudes at resonances 

There was no need to analyze the coupling terms F m in the previous 
section to determine the resonant voltages, but it becomes necessary when 
the AC components A m are wanted. These mode amplitudes are important 
information for possible oscillator applications because one can then estimate 
the power available from the device. 

One way to calculate the amplitudes is to linearize the coupled modal 
equations (||) with respect to A m . One would obtain a linear algebraic system 
with coupling between neighboring modes only. We have tried this approach, 
but the estimated magnitudes became too large to justify the linearization, 
and they did not agree well with the numerically obtained mode amplitudes. 

Therefore, we have used another approximation that can cope with 
larger amplitudes. Since analysis of the modal equations in general appears 
formidable, we focus only on the steps m = 1 and 2, and make several ap- 
proximations. 

First, as we see in the spectra in Fig. ||(a,b), A3 and higher modes have 
much smaller magnitudes on these two steps, so they will be neglected. Then, 

0( s ) = x + (A ie ix + c.c.) + (A 2 e 2ix + c.c.) 
= x + a\ sin(x + 0\) + a 2 sin(2x + 62) 

where 

A m = -(i/2)a m e ie ™. 

The harmonic balance method, which is commonly used for approximating a 
finite amplitude solution, would neglect A 2 as well. Thus, the ansatz ( |l7|) 
can be thought of an extension of the method to a multi-harmonic case. 



Including two harmonics, however, introduces coupling terms between them, 
and makes the problem substantially harder. By the Fourier-Bessel expansion, 
the coupling is expressed as 

oo 

Fo = J -(2n+i)( a i)Jn(a2)sm[n9 2 - (2n + 1)9%] 



n=— oo 



and 

F m = -(*/2) £ J m -(2n + l)(«l)^(«2)e 4 ^ 2+ ( m - 2 "- 1 ) ei ) 
n 

+ (i/2)£ J_ m -( 2 „ + i)(ai)Jn( a 2)e- i(n< ' 2 - (m+2n+1)ei) 

n 

for m/0, and J m are Bessel functions. 

To make analytical progress we simplify the coupling by taking only the 
most dominant terms (for small a's) and neglecting all the others. This results 
in 

F 1 = -(t/2)Jo(ai)J (o2), ^2 = -(i/2)J 1 (ai)J o (a 2 )e i01 . 

The higher order Bessel functions ^(ai), ^(ai), • • • and Ji(a 2 ), ^2(^2), • • ., 
are neglected, which can be justified in the range, say, |ai| < 2 and 1 0,2 1 < 1. 
We substitute them into @ for m = 1, 2 and obtain 

(J (ai)/ ai ) 2 = (i?+rV)/J 2 ( fl2 ), tan 5i = -Vu/8 1 

and 

(Jo(a 2 )/a 2 ) 2 = (df + 4TV)/J?(ai), tan(# 2 - 9 X ) = -2Tuj/8 2 . 

It is convenient to think of u as the control parameter and determine a%, 
a,2, 9\ and 9 2 . The arguments of tangent are defined between and n. As 
ui increases from 0, 5\ and 5 2 decrease monotonically, and changes sign from 
positive to negative at the resonance frequencies. Thus, 9\ decreases from n 
to 0, with the crossover 9\ = w/2 at u = u\. Similarly, {9 2 — 9\) decreases 
from 7r to 0, and crosses n/2 at w = uj 2 /2. 

The amplitude equations of (22,23|) are more complicated. It follows from 



( p3| ) that a 2 = when a\ = 0. That is, the second mode is excited only in the 
presence of the first mode. On the other hand, the first mode can be excited 
alone since a\ > in (|22|) when a 2 = 0. The feedback from the second mode 



through Jo(a 2 ) in ( f22| ) is not essential when considering a\. 

To extract more information we make a further approximation. We first 
linearize the Bessel functions on the right hand sides, i.e. Ji(oi) ~ a%/2 and 
Jo( a 2) ~ 1. (Consequently, Fi becomes independent of A 2 , and the feedback 
from the second mode is neglected. Then, (|22|) reduces to the equations 
obtained through the method of harmonic balance 

On the left hand sides of ( f22^]2"3| ) , the function Jo(a)/a is monotonically 
decreasing from 00 to zero as a increases from to the first zero of Jo, which 
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is a* ~ 2.4. Thus, for any u>, a can be found within and a* . Approximating 
the function crudely by Jo(a)/a « 1/a would make the range of a unbounded. 
This is not desirable near the resonances. We could include the next term 
from the expansion of Jo, but this would make the following algebra more 
complicated. Alternatively, we replace the function by another one which has 
the same asymptotic behavior as a —* and a bounded domain for a. We use 

(J (a)/a) 2 « 1/a 2 - 1/4. 

This has a slightly different saturation amplitude at a* = 2 but overall the 
approximation is good; the error in a introduced by this replacement is less 
than 5% for a < 0.7 and at most 20% when a saturates. In return for the error 
introduced at this stage, we obtain simple expressions for the amplitudes as 

ai = (bj + 1/4)- 1 / 2 and a 2 = (b 2 2 + 1/4)" 1 / 2 , 



where 



b\ = 5\ + TV and b\ = 4($f + 4rV)/a 2 . 
Once the phases and amplitudes of the modes are determined, the instan- 



taneous voltage is obtained from the time-derivative of (17). In our normal- 
ization, 

V/I c R n = Tuj [1 + ai cos(x + 6>i) + 2a 2 cos(2x + 2 )] • 
Thus, the DC voltage and the AC amplitude of the mode m are 
Voc/IcRn = and Vac/ I c Rn = Tuma m . 

In Fig. 0, we compare the amplitude estimate with simulations. On the 
first two steps, the peak AC voltage amplitudes for the first two Fourier modes 
are obtained from simulations, shown as data points. The estimates of the 
resonance frequencies uj = uj\ and w 2 /2 from ( |l3| ) are used. The theoretical 
curves show a reasonable agreement with the data on the Eck step. On 
the second step the comparison is worse, but the overall magnitude and f- 
dependence are estimated fairly, considering the several approximations made 
during the derivation. 

The maximum AC voltage oscillation of the first mode on the Eck step is 
achieved at / = 0.5. The maximum of the second mode on the second step 
happens near the lower critical frustration (/ ~ 0.2). The vertical axes in 
Fig. |7| are shown in terms of AC voltages. Converted back to the amplitudes, 
a\ < 1.75 and a 2 < 0.4 on the Eck step for 0.2 < f < 0.8. This is thought 
to be within the validity of the assumptions. On the second step, a\ < 1.8 
for 0.15 < f < 0.45, but a 2 exceeds unity when / < 0.25. The bending of 
the theoretical curves in (b) when / < 0.25 may be attributed to the large 
predicted magnitude of a 2 . Including terms involving J\(a 2 ) in the modal 
equations would be necessary to make a better estimate in the region. 
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FIG. 7. Amplitudes of the first two AC voltage harmonics on (a) the m = 1 (Eck) step and (b) 
the m = 2 step vs. /. The theory (25) is represented by a solid (dashed) line for the first (second) 
harmonic. Filled circles and open squares show the amplitude of the first and second harmonics, 
respectively, from simulations. The values of Aj and T are the same as in Fig. ||. As expected, the 
amplitude at the second frequency is quite small on the m = 1 step (a), while the amplitudes of 
the first two harmonics are comparable on the m = 2 step (b). 



Generally, our assumption is violated when a\ saturates in (|25|), or when 
b\ becomes small. On the Eck peak b\ ~ (2rAj sin(-7r/)) 2 , so the assumption 
breaks down for very small V and/or Aj. Similarly, a 2 exceeding unity indi- 
cates that we need more coupling terms included from (||). On the m = 2 
step, a 2 < 1 is fulfilled when r 2 c^[l + Y 2 uj\ + 4(w? - w|/4) 2 ] > 3/4. Again, 
small Aj and/or T make it easier to violate the criterion. Thus, the validity of 
our approximations is limited to when TAj is relatively large and resonances 
are not so strong. 



C. Power at resonances 



The calculation of the mode amplitudes provides a theoretical upper limit 
for the power available from the Josephson oscillator at resonance. If we 
assume a matched load condition, then the power dissipated in the load is 
just the power dissipted in the matched resistance. Using the RSJ model and 
assuming sinusoidal voltage oscillations, this correponds to P = V^ c /8R arr , 
where R m -r = Rn/N is the array resistance. The power available from mode 
m is then 

Pm N (V AC \ 2 N 2 

JlR- n = j{7jrJ =Y^ mam) (29) 

with a m estimated in (p5|). Therefore, the vertical axis of Fig. |7] is also 
^j8P m /Nl£R n , and the /-dependence of the powers can be read directly from 
the figure. When biasing on the Eck step, the largest power and the highest 
frequency are obtained at / = 0.5. When the system is biased on the second 
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step, the largest power is obtained near / = 0.2, and the highest frequency is 
achieved near / = 0.25. 

We take our experimental values (I c R n = 0.93 mV and R n = 16.6 0) and 
the maximum amplitudes in Fig.]?] (Vac,i = 0.65/ c i? n on the Eck step and 
Vac, 2 = 0.36I c R n on the second step). This predicts a maximum of about 
150 nW from the first harmonic at the first resonance and 46 nW from the 
second harmonic at the second resonance. These numbers are much smaller 



than previous estimates [24|, and, according to measurements [ 25 , 26 1 , are 
probably more realistic. 

Our simulation and analysis show that multiple Fourier modes are excited 
on the new steps m > 1. For possible applications as oscillators, this is 
not desirable since the AC power is therefore distributed among the modes, 
instead of being concentrated in one mode, as at the Eck step. Usually in such 
a case one could simply increase the drive, to increase the output at the desired 
frequency. However, in this system the driving current is limited because if 
it is increased beyond the top of the step, the resonance becomes unstable. 
Furthermore, since the maximum output power is small for a single row, 
methods are needed to combine power from multiple rows while preserving 
the frequency content of the resonance. 

Both discrete rows of underdamped Josephson junctions and continuous 
long junctions operating at the Eck step have already been proposed as high 
frequency oscillators [1C,27|. Although no single experiment has directly com- 
pared the two systems, studies indicate that the ouput power levels are com- 
parable, while the output impedence of the discrete system may be advanta- 
geously higher In order to increase power levels in either case, additional 
oscillators are required. Stacks of continuous long junctions have been fab- 
ricated and measured [28-30], as well as the analogous inductively coupled 
discrete arrays ||||. However, in both systems, the power is only increased 
if the oscillating elements are phase-locked and in-phase. The existence of in- 
phase as well as anti-phase states have been well established for both discrete 
and continuous systems, through numerics [31|, simulations |2.32], and exper- 
iments j5,33|. Unfortunately, in both cases the anti-phase pattern is preferred 
at the Eck step |5|,^3|, and the output oscillations of adjacent elements almost 
completely cancel. In addition, the difficulty (and necessity) of fabricating 
identical oscillators in the case of continuous stacks limits the possibility for 
improving this system. 

On the other hand, the presence of higher harmonic resonances in discrete 
systems may actually be an advantage in terms of power combining. As 
observed in || , when the discrete system is biased on the stable m = 2 peak, 
the fundamental oscillations of the rows are anti-phase (shifted by 7r), but 
the second harmonics are still in-phase. Thus, the output power from the 
second mode is enhanced. It is important to note that such modes do not 
exist in uniform continuous long junctions. So although the power level of a 
single row biased in this state is lower than the oscillators mentioned above, 
the potential for power-combining in this state is much greater. Simulations 
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indicate that this in-phase resonance for the second mode is stable for many 
coupled rows @, which is promising for oscillator applications. 

V. SUMMARY 

We have observed experimentally and through simulations new resonance 
steps in long one-dimensional Josephson junction arrays in the underdamped 
and discrete regime. They emerge in the I—V characteristic below the Eck 
voltage, and are tunable in the magnetic field /. We have shown that the 
resonance mechanism can be explained by neglecting the boundary effects 
and assuming a traveling wave solution. Fourier components of the solution 
resonate to create a waveform which appears to consist of a kink and radiation, 
phase- locked to each other. We have employed a two-mode extension of the 
harmonic balance method in order to simplify the coupling of the modes. We 
have derived an analytic formula for not only the resonance voltages but also 
the amplitudes of the modes. Finally, we used these results to predict the 
power available from such Josephson oscillators and discuss its dependence 
on the system parameters. 
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